
package wsc13


import collection.mutable.ArrayBuffer
import math.{ abs, sqrt, pow }

import scalation.process._

import scalation.linalgebra.{ MatrixD, VectorD }
import scalation.linalgebra_gen.Vectors.VectorI
import scalation.metamodel.QuadraticFit
import scalation.random.{ Random, Uniform, Variate, Exponential }
import scalation.queueingnet.JacksonNet 
import scalation.stat.StatVector
import scalation.optima.IntegerTabuSearch
import scalation.util.Monitor

trait EmergencyModelParams
{

    val λ = VectorD (6.0 / 60.0,  3.0 / 60.0, 0.0, 0.0,  0.0)
    val μ = Array(VectorD (15.0 / 60.0, 10.0 / 60.0, 3.0 / 60.0, 7.0 / 60.0, 10.0 / 60.0),
                  VectorD (15.0 / 60.0, 10.0 / 60.0, 3.0 / 60.0, 7.0 / 60.0, 10.0 / 60.0))

    
    implicit def vecIA (v: VectorI): Array[Int] = {
        val a = Array.ofDim[Int](v.dim)
        for (i <- 0 until v.dim) a(i) = v(i)
        a
    }

    implicit def vecI2D (v: VectorI): VectorD = {
        val x = new VectorD(v.dim)
        for (i <- 0 until v.dim) x(i) = v(i)
        x
    } // vecI2D

} // EmergencyModelParams

/**
 * Emergency Room Model
 *
 * @param x       the staffing level
 * @param animate animate the model
 */
case class EmergencyModel (x: VectorI, animate: Boolean = false) 
    extends Model ("UCFModel", animate) with EmergencyModelParams 
{

    val (nArrivalsFrontDoor, iArrivalFRV) = (Int.MaxValue, Exponential(1.0/λ(0)))
    val (nArrivalsAmbulance, iArrivalARV) = (Int.MaxValue, Exponential(1.0/λ(1)))

    val μTN = Array(Exponential(1.0/μ(0)(0)), Exponential(1.0/μ(1)(0)))
    val μRN = Array(Exponential(1.0/μ(0)(1)), Exponential(1.0/μ(1)(1)))
    val μMD = Array(Exponential(1.0/μ(0)(2)), Exponential(1.0/μ(1)(2)))
    val μNP = Array(Exponential(1.0/μ(0)(3)), Exponential(1.0/μ(1)(3)))
    val μAC = Array(Exponential(1.0/μ(0)(4)), Exponential(1.0/μ(1)(4)))

    var Δ = Uniform (0, 0)

    var endTime = 100.0 // * 7.0 // 7 days

    // Turn of monitoring
    Monitor.traceOff

    // Set the staffing level
    val nTN = x(0) // Triage Nurses
    val nRN = x(1) // Registered Nurses
    val nMD = x(2) // Medical Doctors
    val nNP = x(3) // Nurse Practitioners
    val nAC = x(4) // Administrative Clerks

    // Sources
    val frontDoor     = new Source ("frontDoor", this, WalkInPatient, nArrivalsFrontDoor, iArrivalFRV, (10., 180.))
    val ambulance     = new Source ("ambulance", this, AmbulancePatient, nArrivalsAmbulance, iArrivalARV, (10., 85.))
 
    // Queues
    val tnQ           = new WaitQueue ("tnQ", (120., 185.))
    val rnQ           = new WaitQueue ("rnQ", (300., 185.))
    val mdQ           = new WaitQueue ("mdQ", (400.,  85.))
    val npQ           = new WaitQueue ("npQ", (400., 285.))
    val acQ           = new WaitQueue ("acQ", (500., 185.))

    // Resources
    val tn            = new Resource ("tn", tnQ, nTN, μTN(0), (170., 180.))
    val rn            = new Resource ("rn", rnQ, nRN, μRN(0), (350., 180.))
    val md            = new Resource ("md", mdQ, nMD, μMD(0), (450.,  80.))
    val np            = new Resource ("np", npQ, nNP, μNP(0), (450., 280.))
    val ac            = new Resource ("ac", acQ, nAC, μAC(0), (550., 180.))

    // Sinks
    val door          = new Sink ("door", (650., 180.))

    // Transports
    val frontDoor2tnQ = new Transport (   "2tnQ", Δ, frontDoor, tnQ )
    val ambulance2rnQ = new Transport (   "2rnQ", Δ, ambulance, rnQ )
    val tn2rnQ        = new Transport ( "tn2rnQ", Δ,        tn, rnQ )
    val rn2mdQ        = new Transport ( "rn2mdQ", Δ,        rn, mdQ )
    val rn2npQ        = new Transport ( "rn2npQ", Δ,        rn, npQ )
    val md2acQ        = new Transport ( "md2acQ", Δ,        md, acQ )
    val np2acQ        = new Transport ( "np2acQ", Δ,        np, acQ )
    val ac2door       = new Transport ("ac2door", Δ,        ac, door)

    // random number generator
    val r = Random () 
    
    // add components to the simulation model
    addComponents (List (
        frontDoor, 
        ambulance, 
        tnQ, tn, 
        rnQ, rn, 
        mdQ, md, 
        npQ, np, 
        acQ, ac, 
        door,
        frontDoor2tnQ, ambulance2rnQ,
        tn2rnQ, rn2mdQ, rn2npQ, md2acQ, np2acQ, ac2door
    )) // addComponents 

    def sumMeanWaitingTimes () = Seq(tnQ, rnQ, mdQ, npQ, acQ).map(r => r.durationStat.mean).reduceLeft(_+_)
    f = sumMeanWaitingTimes

    var high = 0
    var low  = 0

    abstract class Patient (pType: String) extends SimActor ("p-%s".format(pType), this)
    {
        val pclass = if (r.gen <= 0.25) 0 else 1
    } // Patient

    case class AmbulancePatient () extends Patient ("ambulance")
    {
        def act () 
        {
            ambulance2rnQ.move ()
            if (rn.busy) rnQ.waitIn else rnQ.noWait
            rn.utilize(μRN(pclass).gen)
            rn.release
            pclass match {
                case 0 => {
                    low += 1
                    rn2mdQ.move
                    if (md.busy) mdQ.waitIn else mdQ.noWait
                    md.utilize(μMD(pclass).gen)
                    md.release
                    md2acQ.move
                }
                case 1 => {
                    high += 1
                    rn2npQ.move
                    if (np.busy) npQ.waitIn else npQ.noWait
                    np.utilize(μNP(pclass).gen)
                    np.release
                    np2acQ.move
                }
            } // match
            if (ac.busy) acQ.waitIn else acQ.noWait
            ac.utilize(μAC(pclass).gen)
            ac.release
            ac2door.move
            door.leave
        } // act
    } // AmbulancePatient
    
    case class WalkInPatient () extends Patient ("walkin")
    {
        def act () 
        {
            frontDoor2tnQ.move
            if (tn.busy) tnQ.waitIn else tnQ.noWait
            tn.utilize(μTN(pclass).gen)
            tn.release
            tn2rnQ.move
            if (rn.busy) rnQ.waitIn else rnQ.noWait
            rn.utilize(μRN(pclass).gen)
            rn.release
            pclass match {
                case 0 => {
                    low += 1
                    rn2mdQ.move
                    if (md.busy) mdQ.waitIn else mdQ.noWait
                    md.utilize(μMD(pclass).gen)
                    md.release
                    md2acQ.move
                }
                case 1 => {
                    high += 1
                    rn2npQ.move
                    if (np.busy) npQ.waitIn else npQ.noWait
                    np.utilize(μNP(pclass).gen)
                    np.release
                    np2acQ.move
                }
            } // match
            if (ac.busy) acQ.waitIn else acQ.noWait
            ac.utilize(μAC(pclass).gen)
            ac.release
            ac2door.move
            door.leave
        } // act
    } // WalkInPatient

} // EmergencyModel

object EmergencyModelTest extends App
{
    val x = new VectorI(1, 1, 2, 2, 2)
    val m = new EmergencyModel(x, false)
    val results = m.simulate()
    println("done")
} // EmergencyModelTest 

trait EmergencyOpt extends EmergencyModelParams 
{

    val d = 1.0
    val s = VectorD (10, 50, 120, 80, 8)
    val f = VectorD (80, 150)

    def u (x: VectorI): Double = {
        for (i <- 0 until x.dim if x(i) < 1)  return Double.NegativeInfinity
        for (i <- 0 until x.dim if x(i) > 10) return Double.NegativeInfinity
        _u(x)
    } // f

    // this version is used for multiple linear regression
    def _u (x: VectorI): Double = {
        val m = new EmergencyModel(x)       // create model 
        val results = m.simulate()          // simulate and gather stats
        val w = m.sumMeanWaitingTimes()     // sum of average waiting times
        val n = VectorD (m.low, m.high) // patients served of each type
        val c = x dot s                     // operating cost
        val r = f dot n                     // revenue for patient service
        val p = r - c                       // net profit
        p - d * w * n.mean                  // return overall utility
    } // f

} // EmergencyOpt

object EmergencyOptTest extends App with EmergencyOpt 
{

    val x0  = new VectorI(1, 1, 1, 1, 1)
    val its = new IntegerTabuSearch (u, null, 1)
    var (x, result) = its.maximize (x0)
    println ("optimal solution x = " + x + " with an objective value f(x) = " + result)

} // EmergencyOptTest

class EmergencyJacksonModel (x: VectorI) extends EmergencyModelParams 
{
    //                            (DESTINATIONS)
    //                            TN    RN    MD    NP    AC
    val p  = new MatrixD ((5, 5), 0.00, 1.00, 0.00, 0.00, 0.00, // TN
                                  0.00, 0.00, 0.25, 0.75, 0.00, // RN
                                  0.00, 0.00, 0.00, 0.00, 1.00, // MD (SOURCES)
                                  0.00, 0.00, 0.00, 0.00, 1.00, // NP
                                  0.00, 0.00, 0.00, 0.00, 0.00) // AC

    // the jackson queueing network
    val jqn = new JacksonNet (p, λ, μ(0))
    
    def patientsServed: Double = {
        jqn.rho(4) * jqn.k(4) + jqn.nQueue(4)
    } // patientsServed

    def sumOfWaitingTimes: Double = {
        (0 until jqn.m).map( i => jqn.nQueue(i) / jqn.lambda(i) ).sum
    } // sumOfWaitingTimes

} // EmergencyJacksonModel

object EmergencyJacksonTest extends App with EmergencyModelParams 
{
    val x0 = new VectorI (1, 1, 1, 1, 1)
    val qm = new EmergencyJacksonModel (x0)
    qm.jqn.check
    qm.jqn.report
} //EmergencyJacksonTest

trait EmergencyJacksonOpt extends EmergencyModelParams
{
    val d = 1.0
    val s = VectorD (10, 50, 120, 80, 8)
    val f = 100.00

    def u_q (x: VectorI): Double = {
        for (i <- 0 until x.dim if x(i) < 1)  return Double.NegativeInfinity
        for (i <- 0 until x.dim if x(i) > 11) return Double.NegativeInfinity
        val m = new EmergencyJacksonModel (x)  // create queueing model
        val w = m.sumOfWaitingTimes            // sum of the waiting times
        val n = m.patientsServed               // patients served
        val c = x dot s
        val r = f * n
        val p = r - c
        p - d * w * n
    } // u_q

    // this unconstrained version is used for multiple linear regression
    def u_q2 (x: VectorI): Double = {
        val m = new EmergencyJacksonModel (x)  // create queueing model
        val w = m.sumOfWaitingTimes            // sum of the waiting times
        val n = m.patientsServed               // patients served
        val c = x dot s
        val r = f * n
        val p = r - c
        p - d * w * n
    } // u_q2

} // EmergencyJacksonOpt

object EmergencyJacksonOptTest extends App with EmergencyJacksonOpt 
{
    val x0  = new VectorI(1, 1, 1, 1, 1)
    val its = new IntegerTabuSearch (u_q)
    var (x, result) = its.maximize (x0)
    println ("optimal solution x = " + x + " with an objective value f(x) = " + result)
} // EmergencyJacksonOptTest

class EmergencyRegressionModel (u: (VectorD) => Double) extends EmergencyModelParams
{

    // construct a quadratic fit object   
    val qf = new QuadraticFit (u, 5)

    // the current input vector
    var current: VectorI = null

    // the design grid
    var grid = null

    // returns true if x is in the current quadrant
    def sameQuadrant (x: VectorI): Boolean = {
        if (current == null) false
        else (0 until x.dim) map { i => i match {
                case _ if (x(i) <= 5 && current(i) <= 5) => true
                case _ if (x(i) >  5 && current(i) >  5) => true
                case _ => false
            } // match
        } reduceLeft (_ && _)
    } // sameQuadrant

    // returns the design grid for the quadrant containing x
    def designGrid (x: VectorI): ArrayBuffer [VectorD] = {

        //if (curr != null) {
            // TODO only modify the dimension that changed
            //      otherwise we're wasting time
        //} // if

        var grid  = ArrayBuffer.empty[VectorD]
        var parts = Array.ofDim[VectorI](x.dim)

        // each dimension has 3 points depending on where x is
        for (i <- 0 until x.dim) {
            if (x(i) <= 5) parts(i) = new VectorI (1, 3, 5)
            else parts(i) = new VectorI (7, 9, 11) 
        } // for

        // recursively accumulates vec while exhausting rest, adding vec to grid
        // once rest is exhausted 
        // def populate (rest: Array[VectorI] = parts, vec: VectorI = new VectorI (0)) {
        //     if (rest.length == 0) grid += vec
        //     else for (point <- rest(0)) populate (rest.slice(1, rest.length), vec ++ point)
        // } // populate
        // populate ()

        // some temps to help reduce stack space
        var tmp_grid: ArrayBuffer[VectorD] = null
        var tmp_vec = new VectorI (0)

        // populate the grid using an accumulator (similar to the recursion above)
        while (parts.length > 0) {
            if (grid.size == 0) for (point <- parts(0)) grid += tmp_vec ++ point
            else {
                tmp_grid = ArrayBuffer.empty[VectorD]
                for (point <- parts (0); vec <- grid) tmp_grid += vec ++ point
                grid = tmp_grid
            } // if
            parts = parts.slice(1, parts.length)
        } // while

        // return the design grid
        grid

    } // designGrid

    def u_r (x: VectorI): Double = {

        for (i <- 0 until x.dim if x(i) < 1)  return Double.NegativeInfinity
        for (i <- 0 until x.dim if x(i) > 10) return Double.NegativeInfinity

        // if we're in a new quadrant, reform and refit
        if (!sameQuadrant(x)) { 
            // get the 3^5 points to be used for fitting

            // create the design points/vectors at which function f is to 
            // be evaluated
            val _grid = designGrid(x)

            //println("grid vectors")
            //for (vec <- _grid) println (vec)

            // create a design matrix and response vector
            val design = qf.response(_grid)
            
            //println("  design matrix = " + design._1)
            //println("response vector = " + design._2)

            // use multiple regression to fit surface
            qf.fit(design)

        } // if

        // update current coordinate
        current = x

        // use the quadratic regression equation to estimate a value for 
        // the function at x
        val ret = qf.qFormsEval(x.toDouble)
        print("ur(")
        for (i <- 0 until x.dim - 1) print("%d, ".format(x(i)))
        println("%d) = %f".format(x(x.dim - 1), ret))
        ret

    } // u_r

} // EmergencyRegressionModel

/**
// Regression Model for Healthcare Application
class ERRegressionModel (dim: Int, u: (VectorD) => Double)
{   
    val qf = new QuadraticFit (u, dim) // construct a quadratic fit
    var prev: VectorI = null           // last point evaluated

    // returns true if x is in the current quadrant
    def sameQuadrant (x: VectorI): Boolean = {
        if (current == null) false
        else (0 until x.dim) map { i => i match {
                case _ if (x(i) <= 5 && current(i) <= 5) => true
                case _ if (x(i) >  5 && current(i) >  5) => true
                case _ => false
            } // match
        } reduceLeft (_ && _)
    } // sameQuadrant

    // populates a design matrix
    def populate (g: ArrayBuffer [VectorD], v: Array[VectorI], w: VectorI) {
        if (v.length == 0) g += w
        else for (point <- v(0)) populate (v.slice(1, v.length), w ++ point)
    } // populate

    // returns the design grid for the quadrant containing x
    def designGrid (x: VectorI, 
                    lo: Int, hi: Int, 
                    gr: Int, ng: Int): ArrayBuffer [VectorD] = {

        val ngrid = ArrayBuffer.empty[VectorD]
        var parts = Array.ofDim[VectorI](x.dim)

        // each dimension has ng-many points
        for (i <- 0 until x.dim) {
            parts(i) = new VectorI(ng)
            val gsize = lo / gr + hi / gr 
            var group = gsize
            val step  = gsize / ng
            while(x(i) < group) group += gsize
            for (j <- 0 until ng) parts(j) = group + j * step
        } // for

        // populate the design grid
        populate (ngrid, parts, new VectorI (0))

        // return the design grid
        ngrid

    } // designGrid

    def eval(x: VectorI): Double = {
        
        // if we're in a new quadrant, reform and refit
        if (!sameQuadrant(x)) { 
            val grid   = designGrid(x, 1, 11, 2, 3) // build design grid
            val design = qf.response(grid)          // design matrix
            qf.fit(design)                          // fit response vector
        } // if

        prev = x                  // update previous 
        qf.qFormsEval(x.toDouble) // return estimate using regression

    } // eval

} // ERRegressionModel
*/

object EmergencyRegressionTest1 extends App with EmergencyOpt
{
    def f (x: VectorD) = _u (x.toInt)
    val m = new EmergencyRegressionModel(f)
    val x0  = new VectorI (1, 1, 1, 1, 6)
    val its = new IntegerTabuSearch (m.u_r)
    var (x, result) = its.maximize (x0)
    println ("optimal solution x = " + x + " with an objective value f(x) = " + result)
} // EmergencyRegressionTest1

object EmergencyRegressionTest2 extends App with EmergencyJacksonOpt
{
    def f (x: VectorD) = u_q2 (x.toInt)
    val m = new EmergencyRegressionModel(f)
    val x0  = new VectorI (1, 1, 1, 1, 6)
    val its = new IntegerTabuSearch (m.u_r)
    var (x, result) = its.maximize (x0)
    println ("optimal solution x = " + x + " with an objective value f(x) = " + result)
} // EmergencyRegressionTest2


